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We study correlation functions of the Calogero-Sutherland model in the whole range of the inter- 
action parameter. Using the replica method we obtain analytical expressions for the long-distance 
asymptotics of the one-body density matrix in addition to the previously derived asymptotics of the 
pair-distribution function [D.M. Gangardt and A. Kamenev, Nucl. Phys. B, 610, 578 (2001)]. The 
leading analytic and non-analytic terms in the short-distance expansion of the one-body density ma- 
trix are discussed. Numerical results for these correlation functions are obtained using Monte Carlo 
techniques for all distances. The momentum distribution and static structure factor are calculated. 
The potential and kinetic energies are obtained using the Hellmann-Feynman theorem. Perfect 
agreement is found between the analytical expressions and numerical data. These results allow for 
the description of physical regimes of the Calogero-Sutherland model. The zero temperature phase 
. diagram is found to be of a crossover type and includes quasi-condensation, quasi-crystallization 

and quasi-supersolid regimes. 
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I. INTRODUCTION 



There is an ongoing interest in correlation properties of the Calogero-Sutherland model (CSM). From the theoretical 
perspective the Calogero-Sutherland model provides a rare example of exactly solvable model with relatively simple 
structure of eigenstates. In the present work we address the question of coherence properties of the CSM measured 
by the off-diagonal correlation functions. Our study complements the previous results on the diagonal density-density 
correlation properties of the CSM and alows to draw conclusions about the coexistence of both types of long-range 



correlations in an appropriate interval of the interaction parameter. 

Introduced in Ref. p}, the CSM describes a system of particles interacting with a scale-free potential, which is 



inversely proportional to the square of the distance between particles. The ground state wave function was shown to 
have a form which is factorizable over pairs of particles. Each factor is proportional to (xi — Xj) with Xi,Xj being 
the coordinates of particles forming a pair and the parameter A being directly related to the interaction strength. In 
addition to the convenient description of the excitations in terms of non-interacting particles with fractional statistics 
, this particular form of the ground state suggested a possibility to study correlations functions of this model. 
, Indeed, it was noted in pj that for three special values A = 1/2, 1 and 2 the ground state probability of the CSM 
7-H ' coincides with the probability distribution of the eigenvalues of unitary random matrices, so the early results of 
Dyson Q describe the static density correlations in the CSM. The analogy with random matrix theory allows also 
the calculation of dynamical density correlations 3 and dynamical Green's function 0, 0] . To deal with other values 
of interactions the Jack polynomial method has been applied to find correlation functions for integer [7J and rational 
Q values of A. The common drawback of these methods for a rational A = p/q is that the final expression for the 
correlation function is usually given as a sum over fractional excitations involving p + q particle and hole quantum 
numbers, which goes over p + q integrals in the thermodynamic limit. Such a decomposition makes the result appear 
as a highly irregular function of the coupling constant A, leaving little hope of approaching its irrational values. 
Recently, there has been success in obtaining the density-density correlation functions of the CSM for arbitrary 
O ■ coupling A, leading to transparent asymptotic expressions in the long-distance limit. Based on the replica method 
from the theory of disordered systems, this approach involves the representation of the correlation functions using 
a duality transformation as an m-dimensional integral with eventual analytic continuation in to. Similar methods 
were applied recently ^(j to study off-diagonal correlations of one-dimensional impenetrable bosons, equivalent to the 
bosonic CSM for one specific value of A = 1. For recent progress in calculation of correlation functions of integrable 
models such as Heisenberg spin chain see |ll) and references therein. Results for delta-interacting bosons can be found 
in the book [T^. 

In this paper we extend the replica method for studying the off-diagonal correlations (equal-time Green's function) 
and corresponding momentum distributions. The main result is the exact long-distance asymptotic behaviour of the 
one-body density matrix in the form of the Haldane's universal hydrodynamic expansion [13j. We consider both 
bosonic and fermionic statistics of the particles encoded in the symmetry of the wavefunctions using the definitions of 
the original work of Sutherland Q (see also Being irrelevant for the density correlations, the quantum statistics 
affects drastically the results for the off-diagonal correlations. 
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In addition, we study the short-distance behaviour of the one-body density matrix and find several leading terms 
(analytic and non-analytic) in the short-distance expansion. The latter is directly related to the high-momentum tails 
of the momentum distribution. The short-distance (high momentum) physics enters the expressions for potential and 
kinetic energies which we calculate as a function of interactions by using the Hellmann-Fcynman theorem. 

To check our predictions we use a Monte Carlo method to calculate numerically the correlation functions for 
arbitrary distances and different values of the interaction parameter. The advantage of the Calogero-Sutherland 
model is that its ground state wavefunction is known explicitly and can be easily sampled by Metropolis algorithm. 
This permits us to obtain unambiguous results for intermediate distances, where analytical methods fail. 

Combining these results with previous knowledge of the diagonal two-body correlations (pair distribution function) 
allows us to describe different physical regimes of the CSM at zero temperature. In particular, we discuss long- 
and short-range order as a function of the coupling constant A for diagonal and off-diagonal correlations. Since the 
true long-range order is absent in one dimension, we define it through the correlation function of the (local) order 
parameter which has the slowest power-law decay and use the word "quasi" to stress this peculiarity of one dimension. 
We propose a phase diagram which describes (in order of increasing interaction strength) the crossover between three 
physical regimes: the quasi-condensate, quasi-supersolid and quasi-crystal. 

The paper is organized as follows. In Section [H] we introduce the Hamitltonian of the Calogero-Sutherland model, 
its solution for the ground state and define correlation functions of interest. In Section III II we present the calculation 
of the one-body density matrix based on the replica method and discuss the thermodynamic limit in Section IIVI 
Results for the short-distance behaviour of the one-body density matrix and discussion of the kinetic and potential 
energies are found in Section Section IVT1 is devoted to numerical Monte Carlo simulations and discussion of the 
physics of the CSM. In section lVni we propose the phase diagram of the CSM and draw our conclusions. Appendix lAl 
contains mathematical details of our calculations. 



II. THE CALOGERO-SUTHERLAND MODEL 



We consider a finite system of N particles of mass m on a ring of length L. The Hamiltonian of the Calogero- 
Sutherland model is given by the sum of the kinetic energy and pair interactions controlled by parameter A: 



H = 



2m 



N 



(1) 



The interaction between two particles on the ring is inversely proportional to square of the chord distance between 
them. It becomes an inverse square potential in the thermodynamic limit. 

The ground state wave function of the Hamiltonian Q was found in p] and can be written as: 



k<l 



<S>{ Xl ,...,x N ) = C N (\)\[(e 2 * 
Here Cn (A) is the normalization constant given by 

Cft (A) 



ix k /L _ 2-Kixi/L 



i r(i + A) 



JY 



L N T(l + XN) 



(2) 



(3) 



The expression as it stands, is valid only for a particular ordering of particles, for instance x\ < X2 < ■ ■ ■ < %n- 
To extend the expression J5J) to other ordering configuration one has to specify the quantum statistics (bosonic or 
fermionic) of the particles. We modify the expression J2J to take into account the symmetry under permutation of 
the particle coordinates 



= Cn(X)U\ 



J2-nix h /L _ 2-irixi/L 



k<l 



Sb,1 



(4) 



by introducing the factor Sb,f such that for bosons Sb — 1, while for fermions Sf — (— 1) P is the parity of the 
permutation P, where xp 1 < xp 2 < . . . < Xp N . As we shall see different quantum statistics of particles is crucial for 
off-diagonal correlations, which was already noted in the context of the Calogero-Sutherland model in early works pj . 

The main quantity of interest is the one-body density matrix which in terms of the ground state wave function is 
written as 



B.F i 

9i [x- 



y) = N I d N 1 x$ bf (xi,...,x N -i,x)$ b ,f(xi, 
Jo 



,x N -i,y). 



(5) 
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Due to the translational invariance in a homogeneous system the one-body density matrix is a function of the difference 
x — y only. Knowledge of the one-body density matrix enables one to calculate the momentum distribution as the 
Fourier transform: 



n k = j dxe- lkx 9l (x). (6) 

The one-body density matrix has dimensions of density and is normalized so that gi(Q) = n, where n — N/L is the 
particle density. The momentum distribution nk is dimcnsionless, it is defined for ki — 2nl/L, where I is an integer, 
and is normalized to the total number of particles ^2 k n k = Lgi(O) = N. 

In addition to the off-diagonal correlation functions we consider the two-body density matrix (pair distribution 
function). It is defined as 

g 2 (x-y) =N(N-1) J d N ~ 2 x \$(xi, . . . , x N . 2 , x, y)\ 2 . (7) 

The static structure factor is related to g 2 (x): 

S k = l + - I dxe' tkx (g 2 (x)-n 2 ). (8) 
n J 

The two-body density matrix has dimensions of density squared and is normalized so that lim-r^oo g 2 (x) = n 2 , while 
the static structure factor is a dimcnsionless quantity and limfc—xx, S k = 1. It is independent of statistics, as it involves 
only the absolute value of the ground state wavefunction. In the next section we calculate analytically the one-body 
density matrix 101 ■ The results for the pair distribution function Q were obtained in and we reproduce them in 
the Section IVIDl 

III. ONE-BODY DENSITY MATRIX 

To calculate the one-body density matrix we define the dimensionless function Gi(a) such that g\{x) — nGi(2Trx/L) 
and Gi(0) = 1. Factorizing the ground state wave function Q we rewrite the definition JSJ in the form of the average 



27rr(A(l + A/")) 

where the average is defined as 



" = 1 ' N.X 



(f (e'V*, , e«*)>„ iA = fj -j^L | A.(e-)| 2A / (e*,e*, . . . , e^) , (10) 

and Ajv(z) is the Vandermonde determinant 

A N (z) = A(zi, as, ... , z N ) = \\{zi - Zj). (11) 

i<j 

Here we have changed the number of particles from N to A^ + 1 in order to deal with A^ dimensional integrals. This 
difference does not matter in the thermodynamic limit, and for the finite system we will restore the correct number 
of particles in the final expressions. 

To calculate the average in (JSJ we use the replica trick, along the lines of the calculation in |9| and [l0(. Namely, 
consider the following function 

zW(a) = (\{{l~e ie i) m {e ia -e^)™) ■ (12) 

It can be shown along the lines of |l4j that G\{a) is obtained from Zm\a) by the analytical continuation m — > A. 
Later we discuss this procedure in some detail and show how the quantum statistics of the particles, bosonic or 
fermionic, appears naturally in our calculations. For the moment we take advantage of the duality transformation 
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[TH , which enables one to re-express the TV-dimensional integral l|12(l depending on the parameter to as a m-dimensional 
integral depending on N as a parameter: 

-iNma/2 rl m l 

Z£Ht)= Sm{i/X) j( d m x\A m (x)\il[xr\l-x a )i-\l-(l-enx a ) N . (13) 

The duality A «-> 1/A becomes evident by comparing the power of Vandermonde determinants in Eqs. (|10fl . I|12l) and 
Q13[l. Wc put Zm\o) = 1, so the normalization constant is given by the Selberg integral: 

r 1 m i m r 2 (-) v (i + -) 

S m (l/X) = J d™x \A m (x)\i J] xr\l x a )*- 1 = I] r T^Ts r 7^ ■ (14) 



a=l 



The dual representation is an excellent starting point for the asymptotic expansion in the large N limit. In this 
limit the main contribution to the integral over each variable x a comes from the limits of integration x + = 1 and 
x_ = 0, so we expand close to these points: 

x a = X- H — ^ — — o=l, ...A (15) 

N{l-e la ) v ' 

x b = x + — —j- ^ h — — , b = I + 1, . . . ,m. (16) 
To the leading order in N the stationary action is 

(i - (i - e ia ) Xa ) N = ( e r^' _ e ^ = 1 i (17) 

ye e ia , a = I + 1, . . . , to 

and the Vandermonde determinants factorize as 

^ \ 1(1-1) /2 / 1 ^ (m-l)(m-l-l)/2 



AM * U(l-e iQ ) J U(l-e- to ) J A^A^fe), (18) 

which allows the calculation of the remaining fluctuation contributions in the form of the Selberg integrals 0] as 
follows 

z,(i/a) = ia,(oi* ng-v^ = n £^£&+ s . (19) 

Collecting factors arising from the change of variables H16[) and summing over all 2 m saddle points yields the 
following result: 

= E(-l) x( ^ } <( A ) (2°) 

i=o (27V sin f) 

where the factors 

include the combinatorial factor arising from the number of ways to choose I variables x a close to one saddle point, 
x^ say, and to — I variables x in the vicinity of the other saddle point. Having established the asymptotic expression 
(|20|l valid for integer to we consider separately the analytic continuation to — > A for bosons and fermions. 



A. Bosonic statistics 



To obtain the bosonic one-body density matrix one should treat to as an even integer before taking the limit 
to — > A. In this case the main contribution to the sum (|20|) comes from the central point I = to/2. It behaves 
as (iVsin(a;/2))~ m / 2A and substituting to = A yields the result expected from the conformal field theory G\ ~ 
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(iVsin(a/2)) A / 2 . To perform analytic continuation we rearrange the sum by changing the summation index I = 
to/2 + k and letting k run from — oo to +00. The coefficient of the dominant (k = 0) term is given by A^(X) = 

lim m ^A Hm^ 2 ■ The analytic continuation of this expression is described in detail in Appendix [X] The result is 

. / \ \ T 1 /2(i + A) f°° dt t (\ 2(cosh|-l) \ 

The coefficients D\{\) — \im m ^\ H™J 2+k / H™J 2 of the oscillating k ^ terms are obtained straightforwardly: 



k r /W2+a\ k 

D k (X) = lim "TT } \ = TT )} 2 , x} - s (23) 



a=l 



r(^) L\ t (i + j-j 

and can be calculated recursively 



rii+i) n m r(i + i + ^ 



Di(X) = ^,T*' , £> fc+1 (A)= V2 T A k x ' D k (X). (24) 



The correlation function is thus represented as a sum of one smooth component and oscillatory corrections due to 
the short-distance correlations: 

A2r\\ ( 00 cos ( kNa ) \ 

^w-^^+'D-^WteWJ- (25) 

Here we have changed the number of particles back to N, and defined X — (N — 1) sin ^. From the last expression 
one sees that if A is an even integer the infinite sum in the expression l|25|l becomes finite. For A = 2n the correlation 
function contains only n oscillatory components. 



B. Fermionic statistics 



To deal with fermions, m has to be considered as an odd integer prior to the limiting procedure m — * A. In this 
case the dominant contribution to the sum l|2(J|) comes from the terms I = (to ± l)/2, which being combined together 
for m — > A lead to the oscillating behaviour ~ 2 sin (|(iV + /((2iVsin(a/2)) A / 2+1 / 2A of the one-body density 
matrix. This behaviour is again in complete agreement with the conformal field theory. The coefficient of this term 
C 2 (A) = lim m _»> H^r +1 ^ 2 is calculated by an analytic continuation as explained in Appendix 1X1 The result is 

C(X) r !/ 2 (A) r(1/2 + 1/2A) c::p f°° -e"* (---+ — { — * ~ — ^ ) (26) 
° ( j [ ) r(l/2 + A/2) CXP i t £ 4A + (l- e -*)( e */ A -l) )■ [2b) 

and the coefficients F k (\) = lim. m ^\ H^ +k '' 2 / Hm are given by 

_ / A 2 -(2fc-i) 2 y fc - 1 r(i + ^) 
Fk{x) { — a^i — ) llwi + i^y (27) 



or by the recursion relation 



a 2 - (2fc + i) 2 y r(i^ 
A 2 -(2fc-i) 2 J rfr, x 



Fi(A) = 1, ^+i(A) = ( ;' 2 )Z ' ^2 ) ^TTVrM^W- (28) 



Restoring the number of particles and using X = (N — l)sin^, the fermionic one-body density matrix is finally 
given by a sum of oscillatory terms: 

r^iw i °° sin ((k - l/2)Na] 
For A = 2n — 1 the total number of terms in the sum is n. 
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IV. THERMODYNAMIC LIMIT 



We now take the thermodynamic limit of the above expressions by recalling that a = 2ttx/L so that X = 
N sm(irx/L) ~ ttxN/L = nnx = kpx in the limit of infinite N and L. Here we have also defined kp — nn, the 
Fermi momentum by analogy with free one-dimensional spinless fermions. In the thermodynamic limit the expres- 
sions an d EHt become: 



9?{x) 



A 2 {\) 
\2k F x\ x / 2 

C 2 (A) 



2 EM) 



\2k F x\ 



-2T-1 A; 



oo 

E 

m—1 



(-1)' 



, O^(A) cos 2mk F x \ 
\2k F x\ 2m2 / x J 

i-iFmW sin(2m — \)k F x 
{2k F x) 2m ^ m - 1 )/ x ' 



(30) 
(31) 



These expressions agree completely with the results of Haldane based on the universal hydrodynamic theory for 
compressible quantum fluids [l3j ■ The power-law universal decay of correlations with the distance provides the value 
of the corresponding Luttinger parameter K = 1/A. The coefficients A(X) D m {\) (|2*3*|) for bosons and C(X) 

(I26|l . F m (X) i|27[) for fermions, are model-specific non- universal numbers. They are strikingly similar to the analogous 
coefficients which appear in the correlation functions of the Heisenberg spin chain [lTl ] . Here we derive them for the 
first time in the case of the Calogero-Sutherland model. 

From the expressions (|30|l , l|31l) we extract the singular behaviour of the momentum distribution using the definition 
©• For bosons, the leading term in l|30|l yields the leading divergence for small momenta: 



n k ~ — ^(A)r(a)cos — 

7T 2 



k 

k F 



(32) 



where the critical exponent a = a(X) is defined as a = 1 — A/2. This result is valid for < A < 2, so that < a < 1. 
For the special value of A = 1, which corresponds to the system of impenetrable bosons we have a = 1/2. For 
the value A = 1/2 Sutherland [j| was using a numerical estimate to suggest a = l/\/2 = 0.707. . .. We report here 
the exact value a = 3/4. For A — 2 the critical exponent a = is consistent with the logarithmic divergence found in 
0. The oscillating terms in H30|) contribute to the weaker singularities at the points k — ±2k F , ±4k F , . . .. 

For fermions, due to the oscillating character of the one-body density matrix the dominant singularities appear at 
k = ±k F . The leading behaviour is extracted from the first term in l|31|l . For k — > kp we have 



n k -~ 



1 2-' 3 C 2 {X) 



7T/3 



r(l -/3) cos 



tt(/3-1) 



k — k f 



sgn (kp — k) 



(33) 



and similarly for k — > —k F . The critical exponent is defined as (3 = /3(A) = (3(1/X) = A/2 + 1/2A — 1, so that 
the power law behaviour of the momentum distribution is the same for A and 1/A. The result l|33() is valid for 
< < 1 which implies 2 — \/3 < X < 2 + \/3. For A = 1 we have free fermions and [3(1) = corresponding to the 
Fermi-Dirac distribution. In other the special cases we have /3(l/2) = (3(2) — Additional singularities exist at 
k = ±3k F , ±5k F , . . .. 



V. SHORT-DISTANCE CORRELATIONS OF THE CALOGERO-SUTHERLAND MODEL 

The results in the previous section provide the long-distance behaviour of the one-body density matrix for the 
Calogero-Sutherland model and consequently the small-momenta behaviour of its momentum distribution. It is also 
possible to extract the short-distance properties of this model. We use the method introduced by Olshanii and Dunjko 
[l9| to relate the tails of the momentum distribution to singularities of the wave function. To present this method we 
consider bosonic statistics and discuss later the corresponding modifications for fermions. 

The momentum distribution @ is rewritten in the following form: 

n k = n J dx 2 ■ ■ .dx N \$(k,x 2 , ■ ■ ■ ,x N )\ 2 , (34) 
where we have defined the Fourier transform of the ground state wave function @ with respect to its first coordinate: 

$(k,X2,...,Xff) — / dxe~ ikx $(x, X2, ■ ■ ■ , xn) 
Jo 
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cv(A) n 



\ 2irx k /L _ e 2irxi/L 



2<k<l<N 



L 

dxe~ lkx 

1=2 



n 



e 2nx/L _ ^KXijL 



(35) 



The remaining integral can be evaluated using the following property of Fourier transforms. Let f(z) have a singularity 
of the form f(z) = \z — zo\ a g(z) where g(z) is regular at z = zq and a > — 1, 0, 2,4, .... Then 



dze- ipz f(z) = 2 cos 



ir(a + 1) 



F ( 1 + Q ) | p |a+i +° 



(36) 



If /(z) has several singularities, then the right hand side equals the sum of the corresponding contributions. Using 
this fact and (136*1) the large-fe behaviour of the wave function (|3l>|l is determined by the singularities at the positions 
of the remaining particles: 



2ttV ^(A + l) 
l/cl^ 1 V L 

N 

x *^ e -ikXj J~J \^ e 2irxj/L _ e 2-rrxi/L 
3=2 l& 



T(l + \)C N (\) \\ \e 



2-iTx k /L _ ^2-rrxi/L 



2<k<KN 



(37) 



Substituting this expansion into l|34|) and keeping only fc-independent diagonal terms in the double sum leads to the 
following result for the asymptotic behaviour of the momentum distribution: 



2+2A 



n k ^ K N (X) 

with the constant Kn(X) defined by the following average in the ground state of N — 1 particles: 
K N (X) = 



(38) 



2 2+2A o / tt(A + 1)\ r 3 (i + A)T(1 - A + XN) 



■ COS 



N 2X T(1 + XN) 



f 2 " I 

/ dx 2 ■ ■ -dx N |$(x 2 , • ■ • ,x N )\ 2 Y]_ e 



27TX2/L ^277X1 /L 



2 A 



The result (|38|l is valid also for fermions for A ^ 1, 3, . . ., with a modification of the proportionality constant: 

V(A+1)' 



K N (X) = tan 2 



i^jv(A). 



(39) 



(40) 



The result (I38|) suggests the following short range expansion for the one-body density matrix as a sum of analytic 
and non-analytic functions: 



= l + Cl (A)fc F x+^-^(fc F x) 2 + ^^(k F x) 3 + 



5iO) _ , , _ , c 2 (A) 

2! ' 3! 

a(A)|/fc F x| 1+2A + 0(|A :F a ; | 2 + 2A ). 



(41) 



The coefficients q(A) in the Taylor expansion of the analytic part of g\ are the corresponding moments of the 
momentum distribution 



Q(A) 



dk 

n l / 27m 



fc rife. 



(42) 



Due to the time-reversal symmetry the momentum distribution is an even function, gi(x) is real and odd moments 
vanish. The non-analytic part of the one-body density matrix starts as |a;| 2A+1 with the coefficient which can be 
related to the high-momentum tails (|38|l using l|36|) : 



a(A) 



7TK N {\) 



2cos(7r(l + A))r(2 + 2A) 



(43) 



Non-analyticity of the density matrix at x — reflects the fact seen already from (|38|l that the Z-th moment q(A) 
diverges for I > 2 A + 1. 
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Provided they exist, the moments C21 are the same for bosons and fermions, a fact noticed by Sutherland in [jj. 
Therefore for sufficiently large A the fermionic and bosonic one-body density matrices have the same leading Taylor 
expansion. This fact is easily explained on physical grounds, noticing that the strong repulsion prevents all exchange 
effects and particles do not "feel" the quantum statistics. 

The second moment C2(A) can be calculated explicitly since it is given by minus the kinetic energy per particle (in 
units of the Fermi energy). To obtain it one uses the Hellman-Feynman theorem to extract the potential energy from 
the dependence of the total energy on the strength of the particle-particle interaction. The ground state energy of 
the Calogero-Sutherland Hamiltonian is known exactly [lj and is equal to 

E (H) n 2 \ 2 h 2 n 2 

— = - — - = . (44) 

N N 6 m v ; 

We note that the potential energy is linear in g = A(A — 1) and can be obtained by differentiating the ground state 
energy: 

Eyot 9/d \ gd ir 2 \ 2 (\ - I) h 2 n 2 

-W-N\dg- H j- Nd~g { H) ~ 3(2A-1) m' A > 1/2 " (45) 
The kinetic energy is then obtained as a difference of the total ground state energy l|44|l and potential energy l|45|l : 



Ekin = E- E pot _ tt 2 X 2 h 2 n 2 
N N 6(2A-1) m ' 

Finally, for the second moment C2 we obtain the following expression: 



A > 1/2. (46) 



C2 - -3(2X31)' (4?) 

We plot the total, kinetic and potential energies in Fig.GJ The total energy is always positive and the compressibility 
also remains positive so that the system is (thermo-) dynamically stable. The case A > 1 corresponds to repulsion be- 
tween fermions while for A < 1 the interaction between fermions is attractive. The ground state energy is independent 
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of statistics, yet the interpretation for bosons is more involved: the ground state wave function always describes 
repulsion between particles, despite the negative value of the potential energy for A < 1. This paradox is due to the 
singular character of interactions, namely the requirement that wave functions are zero for coinciding positions of the 
particles. 

As can be seen in Fig. both the kinetic and potential energies diverge as A — > 1/2. It corresponds to the critical 
value of the coupling constant g = A(A — 1) = —1/4, below which particles fall to the center 0>|2O|. For A < 1/2 the 
expression (|46f) is not valid, since it predicts negative expectation value for the kinetic energy. The same applies to the 
potential energy Q45[). However the total energy is finite and analytic in the full range A > as follows from Eq. 144fl . 
Approaching the value X — 1/2 from above the divergences of kinetic and potential energies cancel each other. In 
the interval < A < 1/2 one cannot evaluate separately the potential energy using the Hellmann-Feynman theorem 
(|45|l due to the singular character of the ground-state wave function in thi regime. To the best of our knowledge, 
this intriguing behaviour has not been noticed in the literature on the Calogero-Sutherland model despite the simple 
analysis based on use of the Hellmann-Feynman theorem. 

VI. NUMERICAL RESULTS AND DISCUSSION 
A. Monte Carlo method 

Quantum Monte Carlo methods have been successfully applied to the investigation of the equation of state and 
correlation functions in a number of one-dimensional systems pll \'2'j 123) . We resort to the Quantum Monte Carlo 
technique in order to evaluate multidimensional integrals Eqs. Il">l7l) numerically. This calculation of the correlation 
functions is exact in a statistical sense, the small statistical uncertainty can be reduced by increasing the length of the 
simulation runs. An advantage of the CSM is that the ground state wave function is known exactly and can be written 
in a simple explicit way, as given by Eq. |0J, thus facilitating the calculations at zero temperature. Another important 
advantage of the CSM is related to the "sign problem" of fcrmionic Monte Carlo simulations. Here the permutation 
term can be factorized leading to a much simpler and efficient code than in three-dimensional systems, where the 
symmetrization of the fermionic wave function leads to the evaluation of Slater determinants. We use the Metropolis 
[24| algorithm for sampling the square of the wave function and generating a series of states (Markovian chain) having 
the desired probability distribution. The correlation functions are then calculated as averages over the Markovian 
chains. The calculation of the one-body density matrix g\{x) is performed by displacing a certain particle (let us 
choose x\ as an example) by a distance x, so that x\ = x% + x and averaging the ratio $>(x' 1 ,X2 1 ...) /^{x\,X2, ■••) = 
s b, f EL^i | sin(a: / 1 — xi)/ sin(a;i — Xi)\ x where sb = 1 and sp — Yii^i s S n ( x i ~ x i)/ s E n ( x i ~ x i)- The average is 
performed with the probability distribution |$(xi, ...,xn)\ 2 - Similarly, we accumulate the pair distribution function 
by following the Markovian chain generated by the Metropolis algorithm and measuring the interparticle distance. 
The momentum distribution and static structure factor are calculated by means of Fourier transforms as defined by 
Eqs. (l(ilSI) . In what follows we present our results for the correlation functions in different interaction regimes. 

B. One body density matrix 

The results for the one-body density matrix are presented in Figs. [21 for Bose-Einstein and Fermi-Dirac statistics. 
Comparing the numerical results to Eqs. I|30() and (|31|l we find that the asymptotic expansion works extremely well 
even for distances of the order of mean interparticle separation. It is important to note that the series l|31fl and H3UI) 
are asymptotic rather than convergent. The larger the distance x the more terms should be summed, while small 
distances are well described with only a few terms in the sum. 

It is easy to see from Eq. (|30|l that the off-diagonal long-range order is absent as gi(x) always vanishes for large x. 
Still, for bosons and for small values of A the leading term g\{x) ~ A 2 /\2ttx\ x / 2 has a slow power-law decay, which is 
a manifestation of a quasi off-diagonal long-range order or quasi- condensation. This behaviour is shown in Figs.^-Eb 
for small values of A corresponding to weak interaction between particles. The smaller A is, the more pronounced is 
the presence of a quasi-condensate. In this regime the one-body density matrix remains significantly different from 
zero even at distances much larger than the mean interparticle distance. The slowly decaying off-diagonal correlation 
in this regime is well described by the dominant term in the expansion i|3U|l . Oscillating corrections, corresponding 
to m > 1 terms in Eq. I|30(l . decay rapidly on the scale of a few interparticle separations and are unimportant in this 
regime. 

In this regime the off-diagonal correlations in fermionic systems are qualitatively different. For fermions, the one- 
body density matrix drops quickly from unity and oscillates around zero. In the fcrmionic case the main contribution 
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to the long range asymptotics comes from the leading oscillating m = 1 term, Eq. 13111 . 
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FIG. 2: (Color online) One-body density matrix gi(x) in the thermodynamic limit for different values of the interaction 
parameter A. Solid lines: bosons (upper line); fermions (lower line). Long dashed lines: long range expansion for bosons, 
Eq. I3U1 . (upper line); for fermions, Eq. 1311 . (lower line). Short dashed line (c,d,e,f): short-range expansion Eqs. I41I47H . Note 
that the fermionic long-range expansion, Eq. 131H . is exact for A = 1 (c). 



With increasing A the oscillating m > terms in the long-range expansion of the bosonic one-body density matrix, 
Eq. 1301) . become important. Those contributions introduce oscillations corresponding to momenta 2mkp, m = 1, 2, 
where kp = mi is the Fermi momentum. This behaviour can be attributed to the short-range order in the Calogero- 
Sutherland model, which becomes more important as the interaction strength increases. The results for intermediate 
values of interactions are depicted in Fig. [2]>12k- 
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For a bosonic system these short-range correlations are a manifestation of "fermionization" . Indeed, for A = 1 the 
bosonic CSM becomes equivalent to a system of zero-range impenetrable bosons having the density correlations of 
free fermions. The correlation functions involving the phase are, however, drastically different and were the subject 
of the classic works of Refs. [2^,|2(J. Recently there has been a revival of interest in this model [H H EH, E3 due 
to realization of impenetrable bosons in several experiments with cold atoms [28j- For the fermionic CSM at A = 1, 
the one-body density matrix is given by a standard expression for one-dimensional spinless fermions: 

gf^sinfc^ A = l. (48) 
n kpx 

It is interesting to note that the asymptotic expansion (|31|) is exact for A = 1. The results for the one-body density 
matrix at A = 1 are shown in Fig. both for free fermions and impenetrable bosons. 

In the regime A > 1 the system enters the quasi-crystal regime. The off-diagonal decay of the one-body density 
matrix is greatly enhanced. The oscillating terms in Eq. I|30|) become relevant and oscillations in gf{x) start to be 
visible (see Figs. |2i-|2]3) which reflects the appearance of the short range quasi-crystal order. For A = 2 the bosonic 
one-body density matrix is known exactly [lj and is given by the expression 

gf (x) _ Si(2fc F x) ^ 
n 2kpx 

where Si(a:) is the sine integral function. The expression (|30|1 gives gf ' jn — n/4:kFX — cos2kFX / 4k F x 2 which coincides 
with the large x expansion of l|49|) . We note that this behaviour of gf (x) is critical between long-range and short- 
range correlations, in the sense that the integral of g\{x) over space diverges for A < 2 and converges for A > 2. The 
description in terms of a quasi-condensate is applicable only in the long-range regime A < 2. 

For extremely strong interactions (for example, A = 10, see Fig.^l the potential energy dominates the total energy 
and quantum fluctuations are suppressed. In this regime the role of quantum statistics becomes irrelevant. This is 
due to the onset of the quasi-crystalline order, in which particles form a local crystal lattice so that the exchange 
effects of quantum statistics are irrelevant. As a result the one-body density matrices for bosons and fermions are very 
similar (see Fig. U) . One finds that the fermionic one-body density matrix is positive in a large region and vanishes 
otherwise. The short-range expansion l|41|l describe quite well both fermionic and bosonic one-body density matrices. 
The quasi-crystal order is better probed by density correlations described in the subsection IVI Dl 

We have investigated the short-range behaviour of gi(x) numerically and we found a good agreement with the 
expression (I41|) . The results are presented in Fig. EJ This figure shows the validity of the short-range expansion for 
different values of A and allows to find the region of applicability of Eqs. I|41I47|) . We see that the analytic short-range 
expansion holds in a larger range in the case of fermionic statistics as the omitted higher order terms of the short-range 
expansion are smaller in this case. 

We point out that the coefficient ci in the short-range Taylor expansion l|41|l has a minimum for A = 1. This is 
clearly seen from numerical evaluation of the fermionic one-body density matrix presented in Fig. 12^.. Interestingly, 
there are pairs of A (for example, A=3/4 and A=3/2) with the same c-i. The numerical results of the short-range 
behavior of a bosonic one-body density matrix are presented in Fig. [^Jd. The dominant ci term in the short-range 
expansion in this case is clearly seen for A = 1,2,3, while for smaller values of A the nonanalytic correction becomes 
comparable to the analytic contribution in the considered range 0.06 < nx < 0.5. 

We estimate the leading non analytic term and find that it goes as |cc| 1+2A both for bosons and fermions. In the 
regime A < 1/2 the non-analytic part in Eq. I|41|l provides the leading contribution 1 — g\{z)jn cx a(A)|z| 1+2A . The 
point A = 1/2 is very special. Indeed, the kinetic energy l|46|l diverges at this point (see Fig. I44f> . thus the Taylor 
coefficient ci is also divergent. However, this divergence is compensated by the divergence in the non- analytical term 
a(A)|a;| 1+2A , which for A = 1/2 is of the same order. We note that for A = 1 the power of the non-analytical term 
becomes integer again and leads to a cubic correction, as was obtained by Olshanii et al. [19J. We prove numerically 
the presence of the non-analytic term proportional to |a;| 1+2A . The coefficient of this term can in principal be obtained 
from a best fit to the numerical data. We note that while the decay law is the same in bosonic and fermionic systems, 
the coefficient a(A) depends on the statistics and this leads to drastically different behaviour (e.g. A = 1/5 in Figs.[3J|. 
The short-range non-analytical behaviour is related to the high momenta tails of the momentum distribution discussed 
in the next subsection. 



C. Momentum distribution 



The Fourier transformation © relates the one-body density matrix to the momentum distribution. The numerical 
results for this quantity are presented in Fig. 0] for the cases of bosonic and fermionic statistics. 
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FIG. 3: (Color online) Short-range behavior of the bosonic (a) and fermionic (b) one-body density matrix gi(x) in the ther- 
modynamic limit for different values of the interaction parameter A. Symbols: results of the Monte Carlo simulations, lines: 
analytic part of the short-range expansion (A > 1), Eqs. I41I47H . 
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FIG. 4: (Color online) Momentum distributions for bosonic (a) and fermionic (b) CSM for different values of the interaction 
strength. Lines: in descending value at k — 1.5fci?: A = 10; 3; 2; 3/2; 1; 3/4; 2/3. 



In a system of weakly-interacting bosons, the power law decay of the one-body density matrix results in the 
divergence of the momentum distribution for small values of momenta: is proportional to | | 1 ag f }i ows 
from Eq. I|32fl . This infrared divergence for small A is reminiscent of Bose-Einstein condensation, i.e. macroscopic 
occupation of the zero momentum state. The infrared divergence is present in the quasi- condensate regime for A < 2. 
At A — 2 the infrared divergence becomes logarithmic. From (|49|) we have an exact result |l( for the momentum 
distribution: 

n B [\^% \ k \< 2k f , A = 2. (50) 
\ 0, \k\ >2k f y ' 

We note that the logarithmic divergence for A = 2 separates the power-law divergence for A < 2 and a regular 
behaviour for A > 2, so that A = 2 is the critical value beyond which the quasi-condensation disappears. 

The momentum distribution for a fermionic CSM is presented in Fig. 0Ja). In the non- interacting case A = 1 
the momentum distribution of fermions is given by the step function at k = ±kp. For other values of A the steps 
are absent and the momentum distribution has a power-law behaviour close to ±kp described by Eq. This 
is a characteristic of a Luttinger liquid behaviour present in the fermionic CSM. For very large A the interactions 
completely destroy the Fermi surface and the momentum distribution is a decaying featureless function. 

The behaviour of the momentum distribution for large momenta contains information about the physics at small 
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length-scales. For instance, as we see from Eq. (|50f> . the momentum distribution of bosons at A = 2 has a cusp at 
k = ±2kp. This discontinuity is a consequence of the oscillating term in the long-distance asymptotics, Eq. (|3L)|) or 
(|49(l and is a manifestation of the short-range order. Even for A < 2 the presence of the short-range order shows itself 
in weaker singularities of the momentum distribution for even multiples of Uf for bosons and odd multiples of Uf for 
fermions. 

Apart from the special integer values of the interaction parameter A (even for bosons, odd for fermions), the 
momentum distribution has a non-analytic high momenta tails decaying as l/fc 2 ( 1+A ). At the special value A = 1/2 
the momentum distribution becomes broad, with tails decaying as 1/|A:| 3 , which leads to a divergent mean kinetic 
energy oc J k 2 rik dk for A < 1/2. In the case of zero-range impenetrable bosons (A = 1) the Eq. I|32[l yields a 1/k 4 
ultraviolet behavior and has been discussed in [19). 

The short-range order already present the one-body density matrix and the momentum distribution becomes evident 
in density correlations, which we consider in the next section. 



D. Pair distribution function 



The short-range order is best probed by calculating the pair distribution function g%{x) defined in Eq. (JJJ. It 
gives the probability of finding two particle separated by a distance x. This function involves the absolute value of 
the ground state function and therefore is identical for bosons and fermions. The inverse square potential of the 
Calogero-Sutherland model prevents two particles from occupying the same position, so that 172(0) = 0. For large 
interparticle separation the density correlations decouple and we have a general result 52(2) - * n 2 , \x\ — > 00. The pair 
distribution function has been studied in |9j . For distances larger than the mean interparticle separation the following 
result has been obtained 

^ = 1 " ^7 + 2 V J^tn cos(2mk FX ), (51) 



2\{k F x) 2 ^ (2k F x) 2m2 / x 



where 



nr(i + a /A) 

* w - s r H) ii * (f ) r2 (?) ■ < 52 > 

II r(l - a/ A) V ' — 

a=l 

We have calculated the pair distribution function by using the Monte Carlo method and results are presented in 
Fig. [5] One sees that for small values of A the pair distribution function goes smoothly from zero at short distances to 
the bulk constant at large distances. This is characteristic for a weakly interacting one-dimensional Bose gas or liquid 
|2l|. As the regime of free fermions (A = 1) is approached some oscillations become visible in the pair distribution. 
For this value of interactions the expression l|51J) gives the exact result: 

gaCg) , sin2 k FX 

It is interesting to note that the powers characterizing the decay of the oscillating terms in the pair distribution, 
Eq. (|5T1) and in the one-body density matrix of bosons, Eq. (|3T]|> are closely related. Employing the language of 
electronic systems these oscillating terms can be referred to as Friedel oscillations. It is a general feature that these 
oscillations become more pronounced as one moves to stronger interactions, since they arise due to the tendency of 
repelling particles to form a quasi-crystalline order. Friedel oscillations in the Calogero-Sutherland model have been 
discussed in |29| and similar oscillating behaviour of the same origin has also been observed in numerical simulations 
of the strongly interacting limit of one-dimensional bosons with delta-like interactions |2l| • 

As interactions are made stronger the amplitude of oscillations become larger and number of visible oscillations is 
increased (see Fig.[3J). Although the quantum fluctuations prohibit the formation of a "true" crystal, this behaviour 
of the Calogero-Sutherland model can still be described in terms of the local crystalline order as was pointed out by 
Krivnov and Ovchinnikov in 30] . In the next section we show how the quasi-crystal order manifests itself in the static 
structure factor. 
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FIG. 5: (Color online) Pair distribution function for different values of the interaction strength A. Lines: in descending order 
of the value at the first peak: A = 10; 3; 2; 1; 1/2, 1/5. 



E. Static structure factor 



Density correlations characterized by the pair distribution function of the previous subsection are conveniently 
probed by measuring the static structure factor JSJl. It is related to the dynamic structure factor S(k, u>) which 
characterizes the scattering cross-section of inelastic reactions where the scattering probe transfers momentum hk 
and energy fuv to the system. In atomic gases it can be measured directly by the Bragg spectroscopy technique. By 
integrating out the u> dependence one obtains the static structure factor St- The numerical results for this quantity 
in the Calogero-Suthcrland model are summarized in Fig. for different interaction strengths. 

The behaviour of the static structure factor for small momenta can be described by a hydrodynamic approach. The 
Feynman formula [3l| Sk = h 2 k 2 /2mtk relates the static structure factor to the excitation spectrum exhausted by 
one branch ek- This is the case, when k is small and the excitations are phonons tk = Kkc with the speed of sound c. 
This leads to linear behavior for small momenta Sk = h\k\/2mc. The speed of sound is related to the compressibility 
X = mc 2 = ndfj,/dn, where there chemical potential is be found from (|4^|) as fx = dE/dN. Thus we have 

&-J*L. ^0, (54) 

in full agreement with the numerical results as can be seen in Fig. |SJ Large momentum excitations behave as free 
particles ek = h 2 k 2 /2m, so the static structure factor behaves as Sk — > 1 in the limit of large momenta. 

In the regime of weak interactions A — > 0, the static structure factor is a smooth function, which goes monotonically 
from for k = to unity for large values of momenta (see A = l/5;l/2in Fig. 0. This behaviour is very similar to 
that of rarefied weakly interacting gases [21|. The critical point at which a cusp at k = 2k F appears is A = 1. For 
this critical value of A the static structure factor can be written explicitly and takes a very simple form: 

S k = ( 517' jf| < 2k / , A = 1. (55) 



1, \k\ > 2k F 

The linear phononic behaviour continues until |fc| = 2kp, where the asymptotic constant value is reached. This special 
behavior is a result of averaging the dynamic form factor S(k,ui) over all excitation branches (see, for example 32]). 

For A > 1 a peak appears in the static structure factor indicating the onset of the quasi-crystalline order. In this 
regime the correlations between particles on short distance scale (of order several interparticle separations) become 
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FIG. 6: (Color online) Static structure factor for different values of the interaction strength A. Dashed lines: low momenta 
phononic behavior, Eq. 154H . Interaction parameter (in order of decreasing slope at small k): A = 1/5; 1/2; 1; 2; 3; 10. 



stronger. Similar physics was observed in numerical simulations of the metastable gas-like state in a short-range 
attractive potential ( "super- Tonks-Girardeau" system [2^) or with dipole-dipole interactions [23l ]. 

By further increasing the strength of interactions, the peak at k = ±2kF becomes higher (A = 2; 3; 10 in Fig- EJ- 
For extremely strong interactions several peaks can be observed (see, as an example, A = 10 in Fig. EJ) - The absence of 
coherence observed in the off-diagonal correlations in this regime allows us to conclude that for extremely large values 
of A the Calogero-Sutherland model behaves as a lattice of classical particles, the one-dimensional Wigner crystal 
studied in [30| . 

VII. DISCUSSION AND CONCLUSIONS 

The present study of the correlation functions of the Calogero-Sutherland model allows us to make conclusions 
about the dominant order, or rather, quasi-order present in the system as a function of the interaction parameter A. 
We concentrate on the case of bosons and identify the following physical regimes at zero temperature. 

For small values of A the system is weakly interacting and stays in a gas-like (or liquid-like) state. There is a 
substantial degree of phase coherence in the system as follows from the slow decay of the off-diagonal correlations. 
For A < 2 the bosonic momentum distribution has an infrared divergence. This quasi- condensate regime is 
reminiscent of Bose-Einstein condensation, where the k — state is macroscopically occupied. The divergence in 
the momentum distribution disappears completely for A > 2, thus at this interaction value the Calogero-Sutherland 
system crosses over from a quasi-condensate to a non-condensed state for A = 2 (see Fig. 0). 

By increasing the strength of interactions one finds the appearance of strong positional ordering of particles (quasi- 
crystal) as indicated by the large amplitude of slowly decaying oscillations in the density-density correlation function 
gz^x). The critical value for this crossover from a liquid to a quasi-crystal state is estimated as A = 1, where 
a singularity in the static structure Sk appears (see Fig. 0). For stronger interactions the particles form a one- 
dimensional Wigner crystal with a dominant crystalline order and absence of coherence. 

The region 1 < A < 2 is very special as the one-body density matrix exhibits quasi off-diagonal long-range order 
while there is a quasi-crystalline order in the static structure factor. We denote this regime, where those two features 
are simultaneously present, as a quasi-super-solid in analogy with the super-solid state 33]. The word "quasi" is 
necessary while talking about the phase diagrams in a homogeneous one-dimensional system, as there is no true long- 
range order and no true phase transitions can take place in such systems [34] . The change of the regimes are actually 




FIG. 7: (Color online) Phase diagram of the Calogero-Sutherland system. See text for the explanation. 



crossovers. Indeed, the long-range asymptotics of the one-body density matrix ll.'illll include both a slow power-law 
decay term as well as oscillating terms for all values of A. The crossover takes place, when one kind of term becomes 
dominant over the other terms. 

Aside from the question of quasi long-range order related to the behaviour of the correlation functions for large 
distances on the scale of the interparticle spacing, we found interesting phenomena on the scale of the mean distance 
between particles. We found intriguing behaviour in the special region A < 1/2 where both the potential and kinetic 
energies are divergent, while the total energy and the compressibility remain finite. The divergence first occurs for 
the critical value A = 1/2 corresponding to the critical value of interactions beyond which the fall towards the center 
takes place |2(ij ]. 

The presented study of the correlation properties of the Calogero-Sutherland model in different physical regimes 
is of a high fundamental interest since this model provides one of the rare examples of integrable systems for which 
the correlation functions can be calculated exactly. Our description holds in physical regimes ranging from a weakly 
interacting gas to a strongly correlated crystal-like phase. We note that in an arbitrary one-dimensional gapless 
system supporting long-wavelength phonons, the ground state wave function can be approximately written as Q 
[2ij| ■ This general description holds at distances where the hydrodynamic approach is applicable. So, the study of 
the Calogero-Sutherland model is important for understanding the long- wavelength properties of systems with linear 
low-momenta excitation spectra (also known as Luttinger liquids). 

Apart from its theoretical importance, the Calogero-Sutherland model is relevant for several realistic physical 
systems. The classic example is provided by the compressible borders of Fractional Quantum Hall droplets [26j. The 
experiments on vicinal crystal surfaces [^3 provide yet another physical realization of CSM. Recently, the particular 
case A = 1 of bosonic CSM has been realized in the series of experiments [2^ with tightly confined cold atomic gases. 

It is important to stress that regardless of the absence of the true long-range order, the dominant quasi-order can 
reveal itself in mesoscopic confined systems due to their finite size. The study of physical regimes in mesoscopic 
systems confined in harmonic potentials are certainly of interest and will be considered in our future studies. Among 
other open questions one can name the time-dependent correlations as well as correlation functions of the Calogero- 
Sutherland model at finite temperature. 
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APPENDIX A: ANALYTICAL CONTINUATION OF A, 



We show how to perform an analytic continuation of the constant 

T(w + 1) t~t p / m + c \ 

r2(m/2 + l)l = l ^ A J 11 p ( 1 + ?n£±a\ r2 



A 2 (\)- r ( m + 1 ) Y\r( m + C \T\ r ( 1 + f) (A1) 
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to to = A. The first factor in the right hand side of this expression can be continued straightforwardly. Consider the 
logarithm of 

, / A9 r 2 (m/2 + 1)\ /m + c\ v^ 2 / / a\ , / ?77/2 + a\ , /m/2 + a\\ , 

lD \ r (m + IT ) = |> r (— j + g ( lnr i 1 + a) - lnr l 1 + - V- J - 21nr \~x j j (A2) 

We use the following integral representation |38| for the logarithm of Euler's gamma function 

lnr(z) = r | ( ^^^r + (* " ( A3 ) 



to represent each term in the right hand side of (|A2|) . Summing finite geometric and arithmetic series under the 
integral we get 



which under replacement m = A yields the result 

The calculation in the fcrmionic case is similar. We have 



_ r(m + 1) r 2 (iffi) n" i r m r (l + f ) 

U ml A i — \ F , /m+n ~ • l A0 J 

A T 2 (22±1) jj^-t- r (! + 2!=^2») r 2 / ^=1+2^ 



2A 



Taking the logarithm of this expression, using the integral representation (|A3|I and putting m = A yields the result 
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